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d Abstract 

a 

i The coupled-cluster method is applied to the spin-1/2 antiferro- 

magnetic XXZ model on a square lattice by employing an approxima- 
tion which contains two-body long-range correlations and high-order 
O four-body local correlations. Improvement is found for the ground- 

state energy, sublattice magnetization, and the critical anisotropy 
C*~) when comparing with the approximation including the two-body cor- 

J^T relations alone. We also obtain the full excitation spectrum which 

Q\ is in good agreement with the quantum Monte Carlo results and the 

*0 high-order spin-wave theory. 

l> 

2 1 Introduction 

j> The coupled-cluster method (CCM) is one of the most precise microscopic 

formulations of quantum many-body theories [1-9]. There is a large number 
of successful applications of CCM to a wide range of physical and chemical 
systems. In particular, the applications of CCM to spin systems on discrete 
spatial lattices have produced one of the most accurate results [10-20] . Sev- 
eral approximation schemes have been developed for the application of the 
CCM to the spin lattice systems. Two such successful schemes are the so- 
called SUBra scheme in which all correlations of any range for up to n spins 



are retained and the localised LSUBm scheme in which m or fewer adjacent 
spin sites over all distinct locales on the lattice are retained. Other high- 
order localized approximation schemes such as DSUBm [19] and LPSUBm 
[20] have also been employed. Up to now, most recent studies have presented 
results for the high-order calculations mainly based on the LSUBm scheme in 
which the long-range order correlations are ignored [13-20]. In this paper we 
present results for the ground and excitation states for an antiferromagnetic 
square lattice by combining the SUB2 and LSUB4 approximation schemes 
(SUB2+LSUB4). Due to inclusion of the two-body long-range correlations, 
we are able to obtain improved results for the ground-state properties, includ- 
ing the critical value of the anisotropy, as well as the full excitation spectrum 
which is difficult to calculate by using the localised approximation scheme 
alone. 

The spin-1/2 antiferromagnetic XXZ Heisenberg Hamiltonian in terms of 
spin rasing s + and lowering s~ operators is given by, 

h = \ Ek +s J + s t4 + 2As ^]> C 1 ) 

(id) 

where A is the anisotropy and the sum on (i,j) runs over all the near- 
est neighbor pairs once. The isotropic Heisenberg model is given by A = 1. 
Classically, the ground-state of Eq. (1) is ferromagnetic, with all spins aligned 
along z-axis for all lattice when A < — 1; for |A| < 1 it is antiferromagnetic 
for all bipartite lattice with all spins are aligned along some arbitrary direc- 
tion in the xy— plane; for A > 1 it is antiferromagnetic with spins aligned 
along (±) directions of the z-axis. The classical Neel ground state with all 
up-spins on one sublattice and all down-spins on the other is chosen to be 
the model state in our CCM calculation. In this article, as before, we use 
index i to label sites of the down-spin sublattice and index j for the up-spin 
sublattice. It is useful to introduce a transformation for the local spin axes of 
one sublattice. This is achieved by rotating all up-spins by 180° around the 
y— axis and hence every spin of the system points down in the Neel model 
state with s z = —1/2. This transformation is given by for all j-sublattice 
operators, s T = s x ^is y — > —s and s z — > —s z . The Hamiltonian of Eq. (1) 
after the rotation is rewritten as, 

H = ~\ J> + 4 + s-sj + 2Asfsf]. (2) 

(i,j) 



The ket and bra ground states of the CCM are given in terms of correlation 
operators S and S a respectively, 

|*> = e 5 |$>, S = J2^iCl (3) 

i 

(q\ = ($\Se- s , S=l + J2^ C i^ ( 4 ) 

where the model state |<&) is the rotated Neel state as mentioned earlier 
with all the spins pointing down, C\ and Cj are the so-called configurational 
creation and destruction operators respectively with the nominal index I 
labeling the multi-spin raising and lowering operators as, 



1 N/2 

/ , SlCj = J—- y j 2 , Sh,i2-i n ,jl,J2-jn S i 1 S i 2 --- S i n S j 1 S j 2 --- S j n , (s) 



n=l h,i2—in,jl,J2—jn 

N/2 
[III ° 



_^ 1 ' 

2^ SlCl = J^ 2_^ 2^ < S h,i2-i n ,h,J2-jn S i 1 S i 2 --- S i n S j 1 S j 2 --- S j n , (6) 

n=l ii,i 2 ...i n ,ji,j 2 ...j n 

with the ket-and bra-state correlation coefficients Sj and 5/ to be deter- 
mined variationally as shown below. We note that the bra-state (\&| and the 
ket-state |\&) are not manifestly hermitian conjugate to one another. The 
normalization conditions (^|^) = (^l 1 ! 7 } = (^l^) = 1 is satisfied by con- 
struction. The ground-state Schrodinger equation, H\^f) = E g \^/), can now 
be written as, 

H\*) = E B \Q), (7) 

where the similarity-transformed Hamiltonian H can be written in terms of 
a series of nested commutations as, 

H = e~ s He s = H+[H,S] + i[[H, S], S] + ■ ■ ■ . (8) 

The expectation value of an arbitrary operator O can be written as, 

= (#|<9|tt) = (<S>\Se- s Oe s \$) = (5({<S/,<S 7 }). (9) 



The correlation coefficients {Sj,Sj} are determined variationally by the fol- 
lowing equations, 

JL = o => ($|C/e" 5 Fe 5 |$) = 0, (10) 

5Si 

^L = => ($\Se- s [H, a+]e 5 |$) = 0. (11) 

oSi 

In the followings we will consider a specific approximation, namely the SUB2+LSUB4 
scheme as defined earlier, by a similar truncation in both S and S. 

2 Ground-state energy for the SUB2+LSUB4 
approximation scheme 

As mentioned in Introduction, the SUB2 approximation retains two-spin-flip 
configurations of all orders. In the SUB4 scheme, additional 4-spin correla- 
tions are also included. We hence write the SUB4 ket-state operators as, 

N/2 N/2 

i,j h,i2,ji,J2 

where 6y and g^^-j^ are the two-spin-flip and four-spin-flip correlation 
coefficients respectively. The full SUB4 scheme equations were obtained be- 
fore [11], but they are difficult to solve. Here we consider the SUB2+LSUB4 
scheme which retains ten local configurations as shown in Fig. 1, in additional 
to the other two-body high-order coefficients of the SUB2 scheme. 

As described in general by Eq. (10), the SUB4 approximation consists 
of two sets of equations, the two-spin-flip and four-spin-flip equations. The 
two-spin-flip equations are given by, 

(<5>\s; sje- SsuB4 He SsuB4 \<5>) = 0, (13) 

from which we obtain the subset of the SUB2+LSUB4 approximation as, 

J2 [(1 + 2A6i + 2b\ + Gi)5 r , p + 2(A + 26i)6 r 

p 

+ G 2 5 r ,p 3a + G 3 6r tP3b + 2_^ br'+p+p br-r'-p ] = 0, (14) 
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Figure 1: The graphical representation of the ten local configurations in Eqs. (17)- 
(20) for the short-range part of the SUB2+LSUB4 scheme. The flipped spins with 
respect to the Neel state are indicated by the crosses. 

where p is the nearest-neighbor index vector with four possible values for a 
square lattice, p is any one of them, G a with a = 1, 2, 3 are defined as, 



G 1 = 2g2 + 2g b 4 + 4gt + 8gi, G 2 = g b A , G 3 = gt + 2gi, 



(15) 



and p 3 are 2D vectors containing p with P3 a =(3p a: ,0), and p 3b =(2p x , p y ). The 
four-spin-flip equations are similarly given by, 



($|s,r s - sj sj,e- SsuB4 He SsuBi \$) = 0, 



(16) 



from which we obtain the following four coupled equations, 

AAg a 4 - 4A6 2 + Ab x g c 4 + 8b x g a 4 + 8bjb b 3 - Ab a 3 g c 4 - 8b\g c 4 - 8b b 3 gf = 0, (17) 

b b 



5Ag b 4 - Abl - 2Ab 1 b a 3 + 8b l3 \ + b a 3 g\ + 2b a 3 gf + 2h(b a 3 ) 2 + Ab 2 b b 3 



+ 4hb a 3 b b 3 - 6b b 3 gf - 2%g d 4 - 2b\g{ - 2b\g\ - b%g\ = 0, (18) 
5Ag c 4 - Abl + 2b l ~ 2AM3 + tt>\b\ + 46x(^) 2 + b x g% - b%g a 4 

+ 8b l9 l + 2b l9 d 4 - b a 3 g d 4 - 3b b g d 4 - b\g d 4 - b\g\ - 2b%g\ = 0, (19) 
5Ag d - Abl ~ 2A6i6^ + bf + b 2 x b a 3 + 36 2 ^ + Ab^f + hb b 3 b a 3 

+ hg c A + 8b l9 d 4 - b a 3 g d A - b b 3 g d 4 - 2b\g{ - \{b%gl + 3b b g c 4 

+ b b 3 g\ + b a 3 g b + %gl + b\g\ + b\g\ + b b g b ) = 0. (20) 

These nonlinear equations for the SUB2+LSUB4 scheme are solved firstly by 
Fourier transformation of Eq. (14) and then by iteration method for Eqs. (17)- 
(20). In particular, Eq. (14) becomes after Fourier transformation, 



7 (q)r 2 (q) - 2KT(q) + G 2l3a (q) + G 3l3b (q) + (G 1 + 26 2 + 2A&! + l) 7 (q) = 0, 

(21) 
which is easily solved with the physical solution, 

r(q) = ^y[l-£(q)], (22) 

where the constant K, and the function E(q) are given by respectively, 

K = A + 2&i, (23) 



E(q) = sjl- A;Mq) - *2W<j)7(q) - A^faMq), (24) 
and where 7(q), 73a (q) and 736 (q) are defined respectively by, 

7(q) = ^(cosq x + cosq y ) } (25) 

7sa (q) = -(cosSg^ + 1), (26) 

73fe(q) = -(cos 2ft. + cos ft), (27) 



with the constants kf, k\ and k\ defined by, 



k\ 



1 + 2A&! + 2b\ + d 
(A + 260 2 



t.2_ 9l 

ft 2 - ^ A , OZ. \2' 



(A + 26^ 



^"l 



G a 



(A + 260 2 ' 



(28) 



In any approximation scheme of CCM, the ground-state energy for the Hamil- 
tonian of Eq. (2) is always given by [11], 



E g = (^\H\^>) 



^V(2&! + A), 



(29) 



where z is the coordination number. In Fig. 2 and Table 1, we present 
numerical results for the ground-state energy as a function of the anisotropy 
parameter A in our SUB2+LSUB4 scheme, together with those of the SUB2, 
SUB2+(yf4 and LSUB4 schemes obtained earlier [11] for comparison. As can 
be seen, the SUB2+LSUB4 results are lower than any of the other schemes. 
Furthermore, the critical value of the anisotropy A c =0.847 beyond which 
the solution of Eq. (22) becomes imaginary, is also improved and closer to 
the expected value of 1 than that of the SUB2 scheme (0.798) or that of the 
SUB2+t74 scheme (0.818). In the high-order LSUBm scheme [16], the critical 
values are obtained as A c = 0.763 and 0.843 for m = 6 and 8 respectively, 
and A c = 1 after extrapolation of m = oo is made. The corresponding value 
of A c in the localized schemes are 0.637 in DSUB10 [19] and 0.766 in LPSUB5 
[20]. The physics of this critical point was discussed in details in Ref. [11]. 



Table 1: The ground-state energy per spin for the 2D spin-1/2 XXZ model in the 
SUB2+LSUB4 scheme for some values of A, together with that of the full SUB2, 
SUB2+gf, and LSUB4 schemes [11]. 
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Figure 2: The ground state-energy per spin as a function of A for spin-1/2 XXZ 
model in the full SUB2, SUB2+g£ and SUB2+LSUB4 schemes. The critical ter- 
minating points for each scheme are also indicated. 



3 Staggered Magnetization 

The staggered magnetization for a general spin quantum number s can be 
defined as, 

1 - 

i 

where / runs over all the lattice sites for our rotated Hamiltonian of Eq. (2). 
In the SUB2+LSUB4 scheme we obtain, 

M = 1 - 2 J2 Mr - 2(M + g\g\ + m + gfgt, (31) 

r 

where two-body and four-body bra-state coefficients b r and g^ are determined 
by the second variational Eqs. (11). We solve these equations for the bra- 
state in similar fashion as for the ket-state, namely by Fourier transformation 
for the two-body coefficients and by iteration methods for the four-body 
coefficients. We leave the details to Appendix and show the results in Fig. 3. 
We find that at the critical A c , M c = 0.649 in our SUB2+LSUB4 scheme, 
compared with M c = 0.663 in the SUB2+g2 scheme and M c = 0.682 in the 
SUB2 obtained earlier [11]. Our SUB2+LSUB4 result is in good agreement 
with M = 0.6138 of the 3rd-order spin-wave results [21], M = 0.614 of the 
series expansion calculations [22], M = 0.615 of the quantum Monte Carlo 
calculations [23] at A c = 1. The highe-order LSUBm scheme with m=8 
produces M = 0.705 at A = 1 before extrapolation and M = 0.616 after an 
extrapolation has been made [16]. The corresponding values of M at A = 1 
are 0.712 in DSUB11 scheme [19] and 0.708 in LPSUB6 scheme [20]. 
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Figure 3: The staggered magnetization for the 2D spin-1/2 XXZ model for the 
full SUB2, SUB2+ 9 f and SUB2+LSUB4 schemes. 

4 Spin-wave excitation spectra 

The excited state in CCM is given by applying an excitation operator X e to 
the ket-state wave function, 



l^ e ) = x e |^) = xv|$), 



(32) 



where X e in general is written in terms of the configurational creation oper- 
ators Cf only as, 

X e = J2x}Cl (33) 

i 

with the excitation coefficient \\. From the Schrodinger equation H\^ e ) = 
E e \^ e ), it is straightforward to derive the following equation for the excitation 
coefficient, 

e.x^mC.e-'^Xy 3 ^), (34) 



10 



where e e = E e — E g is the excitation energy. Here, we consider the spin-wave 
excitations by including only a single spin-flip operator, C\ ~ sf, similar 
to the SUB2 scheme as before [11]. After Fourier transform we obtain the 
energy spectrum in this linear approximation as, 

e e = e(q) = ^zKE(q), (35) 

where K and E(q) are as defined before in Eqs. (23) and (24), and z is the 
coordination number. We present the excitation gap, e(q) at q = 0, as a 
function of A in Fig. 4. As can be seen from the figure, the energy gap in 
the SUB2+LSUB4 scheme is smaller than that of the SUB2 and SUB2+^ 
schemes, implying that the energy gap is reduced in the higher-order ap- 
proximations. For all these three schemes, the energy gap disappears at 
their corresponding critical anisotropy A c . It is interesting to compare our 
results for the energy gap with that of the high-order LSUBm scheme [16]. 
At A = 1 our SUB2+LSUB4 gap value is e(0) = 1.05 while the LSUB4 and 
LSUB8 values are much lower at e(0) = 0.851 and 0.473 respectively. By em- 
ploying an extrapolation, the LSUBm scheme produces an energy gap close 
to zero, corresponding to the SUB2+LSUB4 result at the critical A c . The 
much lower energy gap values away from the critical region by the higher- 
order LSUBm scheme are clearly due to the inclusion of the higher-order 
correlations in the excitation operators whereas we only include the linear 
excitation operators in our calculations as given by Eq. (33) with C] ~ sf. 
However, our SUB2+LSUB4 scheme has an advantage of capable of pro- 
ducing the full energy spectra due to inclusion of the long-range two-body 
correlations as discussed below. 

In Fig. 5, we present our SUB2+LSUB4 results for the spin- wave energy 
spectrum of Eq. (35) at A c together with that of the SUB2 results [1 1] , and at 
A = 1, the results of the linear spin- wave theory (LSWT), the series expan- 
sion calculations (SE) [24], and quantum Monte Carlo calculations (QMC) 
[25]. The spin- wave velocity correction factor to the linear spin- wave theory 
in our SUB2+LSUB4 scheme is given by K c = 1.23, in good agreement with 
1.18±0.02 from the series expansion and 1.21 ±0.03 from the quantum Monte 
Carlo calculations. 



11 
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Figure 4: The excitation energy gap e(0) for the 2D spin-1/2 XXZ Heisenberg 
model as a function of A, for the full SUB2, SUB2+ 5 f and SUB2+LSUB4 schemes. 
The two gap values at A = 1 are given by the LSUB4 scheme (•) and LSUB8 
scheme (o) of Ref. [16] where the high-order excitation correlations are included 
as discussed in the text. 
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Figure 5: The spin-wave excitation spectra for the 2D spin-1/2 XXZ Heisenberg 
model at A c for the CCM (SUB2 and SUB2+LSUB4) results, and at A = 1 for 
the linear spin- wave theory (LSWT), the series expansion (SE) [24], and quantum 
Monte Carlo calculations [25] . The energy spectra in (a) are for q x = q y and those 
in (b) are for q y = 0. 
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5 Summary and Conclusion 

In summary, we have obtained here numerical results for the ground-state en- 
ergy, sublattice magnetization, and excitation energy for the spin-half square- 
lattice antiferromagnetic XXZ model using the SUB2+ LSUB4 scheme of 
CCM. We find that our results for the ground-state properties in general are 
improved when compared with those obtained by the SUB2 or LSUB4 scheme 
alone. In particular, due to inclusion of the two-body long-range-order corre- 
lations, the SUB2+LSUB4 scheme is capable of producing improved results 
around the critical regions of the anisotropy, the excitation gaps at q = 0, and 
the full spin-wave energy spectra. Good agreement for the spin-wave spectra 
is found with the high-order series expansion and the quantum Monte Carlo 
calculations. This is contrast to the recent state-of-the-art calculations of the 
LSUBm scheme using computer algebra, where good results of the critical 
properties have been obtained after an extrapolation in the limit m — > oo is 
made [16-20]. Away from the critical points, the long-range correlations are 
less important and the high-order LSUBm clearly provides better numeri- 
cal results due to inclusion of the high-order local correlations. We believe 
that the different approximation schemes in CCM complement each other 
for a more complete description of the physics of the spin-lattice Hamilto- 
nian model, and in particular the SUB2+LSUBm scheme as presented here 
has the advantage of producing the full excitation energy spectrum. Further 
improvement for the excitation energies away from the critical points can be 
obtained by including the higher-order local correlations in the excitations 
operator X e as demonstrated in the LSUBm scheme of Ref. [16]. It will be 
interesting to apply our SUB2+LSUBm scheme to other models such as the 
spin-1/2 XY model. 
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Appendix 

The ground bra-state in the SUB2+LSUB4 

scheme 

Similar to the ket-state equations, the bra state in the SUB2+LSUB4 scheme 
retains the two- and four-body bra-state correlation coefficients defined as b r , 
and gl, g\, g\ and gf respectively. From Eq. (11), there also are two sets of 
equation for the bra-state coefficients. The first set is obtained by taking 
the partial derivatives of the Hamiltonian expectation H with respect to b r , 
thus, 



OH 

db r 



J2 [(! + fl i + 2 ( A + 2b i)bi ~ 4 J2 i>r'br')6 r , f 

p r' 

j- a 3 o rM%n + a> 3 o rM , ih + d^o rt p 5a + a 5 dr tP5b 



l 3 u r,p 3a 
+ a 5°r,p 5c 



v 3 u r,p 3b 
2(A + 2&i)6 r 



2^2b r ,b r . 



(36) 



where the constants, oi, a 3 , a|, ag, a\ and a\ are given as, 



a x =g\ 



2A6i - 2Ab a 3 + 2{blY + 85i6° + Ab a 3 b b 3 + Sg b 4 ) + ^(-8A6j 

+ 16&x^ + 8g a A + Agl) + ^(-2A6 X + Qb\ - 2Ab b 3 + 86i&| 

+ 4(b\) 2 + g\ + Sg c 4 + 2gf) + sfJ(-2A6i + 36? + 26x6° 

- 2A6* + 6&i&§ + b a 3 b b 3 + A(b h 3 ) 2 + s 4 c + Sgf), 
1 



[ < r l</l + ^gf(2bi + 2b 1 b 



A 



9t 



2gl) - Ual + 92) 



+ gl(-2Abt + 4hb a 3 + Abxb b 3 + g b 4 + 2g 



A)t 



(37) 
(38) 



9t(8b( - 8gl - 8gf) + g b 4 (Abf + Ab x b a 3 - Qg d 4 ) + gt(-2Ah + Ab{ 



+ 8b x b b 3 - 3g%) + -gH-AAbi + 6bf + 2b 1 b a 3 + I6M3 



A 



3flS 



2<?f), 



(39) 
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< = -~g\g\, (40) 

a\ = - \ ~gt (g b A + gt) - ~g\(2g\ + 2g d A ) - g%gi (41) 

a c 5 = - \ ^ {g\ + gt + Ag d ) - ~g\{2g\ + gf) - 2~g\g d A . (42) 

and where the 2D vectors p 5a = (5p x , 0), p 5b = (4p x , p y ) and p 5c = (3p x , 2p y ) 
with the nearest-neighbor vector index p = (p x ,p y ). 

The second set of equations for the bra-state are obtained by taking the 
partial derivatives for H with respect to the four-body ket-state coefficients, 
hence, 

' - = 2b 1+ $(4A + 860 + ~9l(bi - b a 3 ) = 0, (43) 



46i + b\ + gl{-2bl + 5A + 86x) + #t(46i - 46^ - 86^) 



— = 2b 1 + b a 3 + g\{-bt - 2b\ + 5A + 86: + 6f ) - - g{ {b\ 

+ 6^ + 6^ + 6^=0, (44) 
dH 
dg\ 



+ I^(26 1 -6»-6^-6 3 5 -36 3 = O, (45) 



dH 

~~4 



8h + 2b\ + g\{-2b\ - 2b% + 2b% - 66*) + g%-b\ - b\ 
+ 26i - bl - 36*) + g d A {-2bl + 5A + 86 x -b%- 6*) - 8^63 = 0. (46) 



Similar to the solution of the ket-state coefficients, in order to find the bra- 
state correlation coefficients, we obtain Fourier transformation of Eq. (36) 
which is solved together with Eqs. (43)-(46) self-consistently. We rewrite 
Eq. (36) in the following simpler form as, 

J2 [(1 + 01 + 2^6i - 4S)5 r , p + a a 3 5 riP3a + a b 3 6 r>P3b + a a h 5 r , P5a 

p 

+ ald r , P5b + a c 5 8 rtP5c - 2Kb r + 2 J^ hK-r'-p = 0, (47) 

r' 

where K is again defined in Eq. (23) and the constant S is given by, 

16 



Z = ^2b r ,b r ,. (48) 

After Fourier transformation, Eq. (47) reduces to 

(1 + 01 + 2Kk - 4E) 7 (q) + A(q) - 2AT (q) + 2 T (q)f (q)r(q) = 0, (49) 

where T(q) and T(q) are the Fourier transformations of the ket- and bra-state 
coefficients respectively, and the function A(q) is given by, 

A{q) =as73a(q) + a37 3b (q) + a^ 5a (q) + a b 5 j 5b (q) + a c 5 j 5c (q), 

with 73a(q) and 73&(q) as given before in Eqs. (26) and (27) and new functions 
defined as, 

75a(q) = "(COS 5^ + 1), 

75?>(q) = -{cosAq x + cos^), 

75c(q) = -(cos3g x + cos2g y ). 

Using the solution for T(q) of Eq. (22) with the definition for E(q) in Eq. (24), 
the physical solution of Eq. (49) for the bra-state is, 

~ 7J 7 (q) + 2A(q) 

r(q) = AKE{q) ' (50) 

where the constant D is defined as, 

D = 2(l + a l + 2Kb l -A~). (51) 

The value of D can be determined self-consistently as follows. We first rewrite 
Eq. (48) as an integral in Fourier space as, 



r i r 2A(q)i r i i , 

/ - D + —^- — - - - 1 dq. 52 

7 4 L 7 (q) J lE(q) J 4 l ' 



i r i 

7T 



The bra-state coefficient b r is obtained by inverse Fourier transformation of 

f(q), 

g _ 1 r c - ^7(q)+2A(q) (53) 

r '~W e 4^£(q) dq ' (5dj 
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and in particular, b\ is given by, 

1 r J D 7 2 (q) + 2A(q) 7 (q^ 



61 



71" 



dq. 



(54) 



, 4^£(q) 

Combining Eqs. (51), (52) and (54), we obtain the following expression for D, 



D 



-1 



n 



-J 1 -?™ 2 * 



E(q) 



(55) 



where the constant c is given by, 



c = I + 01 + 1, 



with the integral J defined as, 



7T- 



A(q) 7 (q) - 2A(q)/ 7 (q) , 2A(q) 



£(q) 



7 (q) 



<iq. 



(56) 



Using the above self-consistency equations for fo 1; 63, 6|, D and S and by 
iteration method, we obtain the numerical values for g%, g\, g c A and gf of 
the four-body bra-state coefficients. The staggered magnetization is then 
calculated by using Eq. (31). 
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